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ABSTRACT 

We investigate the possibility of constraining the ionization state of the Intergalac- 
tic Medium (IGM) close to the end of reionization (z « 6) by measuring the size 
of the H II regions in high-z quasars spectra. We perform a combination of multi- 
phase smoothed particle hydrodynamics (SPH) and three-dimensional (3D) radiative 
transfer (RT) simulations to reliably predict the properties of typical high-z quasar 
H II regions, embedded in a partly neutral IGM (xhi — 0.1). In this work we assume 
a fixed configuration for the quasar lifetime and luminosity, i.e. tq = 10^ yr and 

= 5.2 X lO^^s^^. From the analysis of mock spectra along lines of sight through 
the simulated QSO environment we find that the H ii region size derived from quasar 
spectra is on average 30 per cent smaller than the physical one. Additional maximum 
likelihood analysis shows that this offset induces an overestimate of the neutral hy- 
drogen fraction, a;Hi, by a factor w 3. By applying the same statistical method to a 
sample of observed QSOs our study favors a mostly ionized {xm < 0.06) universe at 
z = 6.1. 

Key words: cosmology: theory - radiative transfer - methods: numerical - intcrgalac- 
tic medium- cosmology: large scale structure of Universe 
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1 INTRODUCTION 

Studies of the Gunn-Peterson (Gunn & Peterson 1965, GP) 
test have now firmly established that at z > 6 the mean 
volume (mass) weighted neutral hydrogen fraction, a;Hi, is 
higher than 10~^ (10"'') (e.g. Fan etal. 2006a); a more pre- 
cise determination is hampered by the strong sensitivity of 
Lya photons resonant scattering to even tiny amounts of 
H I . Nevertheless, assessing the exact value of xm at z ~ 6 
would represent a crucial information to discriminate be- 
tween different reionization scenarios and identify the nature 
and distribution of ionizing sources. 

The size of H ii regions around high-z luminous quasars 
prior to complete reionization is strongly dependent on xui- 
Previous studies {e.g. Wyithe & Loeb 2004; Wyithe etal. 
2005) have tried to use this observable to improve upon the 
above GP constraints. 

In a simple (pre-overlap) reionization picture in which 
isolated quasar H ii regions expand into an intergalactic 
medium (IGM) with mean neutral hydrogen fraction xui, 
the (physical) size of the H ii region is 



Rd 



1/3 



(1) 



where A'"-, 



are the ionizing photons emission rate and 



the quasar lifetime, respectively, and jih is the hydrogen 
number density. This expression does not account for the 
recombinations occurring in the ionized gas, which can be 
neglected as typical quasar lifetimes {tQ = 10^~* yr) are 
much smaller than the IGM recombination time scales at 
these epochs. Also note that eq. [T]does not hold during the 
relativistic phase of the H ii region expansion, which lasts as 
long as Njt > (47r/3)nHa::Hi(ct)^. For typical high-z quasars 
luminosity this phase could last as much as tg . However (see 
e.g. White etal. 2003), the size of the observed H ii region 
along the line of sight (LOS) to the quasar turns out to have 
exactly the same time evolution as deduced by assuming an 
infinite speed of light. In summary, eq.[T]can be safely used 
to describe the size of the quasar H ii region observed in 
quasar absorption spectra and to put constraints on xm. 

In practice, this procedure is affected by a number of 
possible sources of error, which range from uncertainties in 
the quasar properties (lifetime, luminosity, highly biased en- 
vironment) to intrinsic problems in the practical definition 
of the H II region size from the observed spectra. Despite 
all these problems, the available set of quasar spectra with 
a transmission region inside the GP trough have been ex- 
ploited to constrain xm- Wyithe & Loeb (2004) and Wyithe 
etal. (2005) performed a statistical comparison between ob- 
served and predicted radii, and derived a mean a^ni > 0.1 at 
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z « 6. An independent study by Mesinger & Haiman (2004, 
hereafter MH04) , not involving a direct measurement of the 
H II region size, found arm ^ 0.2 by modelling observed 
properties of the Lya and Ly/3 regions of the spectrum of 
the z = 6.28 quasar SDSS J1030+0524. These results seem 
to indicate a largely neutral IGM at those epochs, somewhat 
at odd with recent interpretations of the 3-yr WMAP data 
(Choudhury & Ferrara 2006, Gnedin & Fan 2006). Indepen- 
dent analysis of the same QSOs spectra by Yu & Lu (2005) 
and Fan etal. (2006a) find much lower values for a^ni. Fur- 
thermore a recent theoretical analysis by Bolton & Haenhelt 
(2006) argues that current observed spectra are consistent 
with a very broad range of xm values. 

In this work we move a step further: by performing a de- 
tailed statistical analysis of mock quasar spectra extracted 
from the combined SPH and 3D RT simulations, we quan- 
tify the confidence level at which xhi can be constrained by 
inverting eq[T] 



2 SIMULATIONS AND RESULTS 

We have performed a combination of multiphase SPH and 
3D RT simulations, in order to predict reliably the geometri- 
cal shape of the H ii region around a typical quasar observed 
at z > 6. High- 2 luminous QSOs reside in rare overdense re- 
gions where the IGM physical properties are highly biased 
(Yu & Lu 2005). This bias has been taken into account by 
using a snapshot centered at zq =6.1 of the G5 simulation 
described in Springel & Hernquist (2003). With its large 
computational volume (100/i~^ comoving Mpc on a side) 
and a particle resolution of 2 x 324'^, G5 allows us to prop- 
erly follow the quasar H ii region volume at a sufficiently 
high resolution. The density field is centered on the most 
massive halo, M^aio ~ 2.9 x W^^Mq. This mass is consis- 
tent with that expected for halos hosting high-z luminous 
quasars (Wyithe & Loeb, 2004). 

The SPH density field has been mapped on a Cartesian 
grid with 128^ cells, in order to perform full 3D RT sim- 
ulations. We used the RT code CRASH (MaseUi etal. 2003; 
Ciardi etal. 2001), which follows the time evolution of gas 
ionization state and temperature. Or spatial resolution does 
not allow to resolve the IGM dumpiness at scales above 
« 0.76 Mpc comoving; this prevent us to correctly account 
for nonlinear clumps like halos and for their effects on the 
overall recombination rate and on shielding. 

The quasar is embedded in the most massive halo and 
we have adopted a Telfer template (Telfer etal. 2002) in the 
energy range 13.6el/ — 42eV as its UV spectrum. For com- 
putational economy He physics was not included in the sim- 
ulations. This might lead to a underestimate of the ionized 
gas temperature and could affect the recombination rate. 
A higher temperature might increase the inner (resonant) 
Lya opacity, essentially due to thermal broadening of the 
line. In order to quantify this effect, we have re-calculated 
some examples of mock spectra along the same l.o.s. used in 
the analysis presented in the paper, increasing the temper- 
ature by 30 per cent 0. Since recombinations are negligible 

^ In a previous work (Maselli etal. 2003) we found that includ- 
ing He typically increases the temperature by 30 per cent for a 
mean density of 1 cm~^. As in this work the density are much 
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Figure 1. Map (upper panel) and cut (lower panel) of the sim- 
ulated xjji (logarithmic scale) across the quasar (located at the 
center of the box). The quasar H II region is clearly identified. 

and we have checked that the corresponding thermal broad- 
ening is not affecting the Lya resonant opacity in a sensible 
way. He inclusion would not change the results shown here. 

The quasar radiation is sampled by emitting Np = 10® 
photon packets. We have tested numerical convergence by 
running lower resolution simulations with 0.57Vp and 0.25A'p. 
Furthermore we assume tg = 10'^ yr and = 5.2 x 10'^'^ 
s~^; the latter values have been adopted after MH04 analy- 
sis of SDSS J1030-I-0524. It is fair to note, though, that the 
estimates of tg and should be considered as tentative 
because of: (i) the degeneracy among A'^-,, tg, xhi and the 
adopted spectral template, (ii) uncertainties in the quasar 
redshift, (iii) possible lensing effects, (iv) IGM dumpiness 
(see e.g. White etal. 2003 for a discussion). However, this 
difficulties have little impact on our study as our main aim 
is to estimate the confidence level at which xm can be ex- 
tracted from obserations. 

Initially, the IGM is in photoionization equilibrium with 
an uniform ionizing background (produced by sources other 
than the considered quasar) with a mean photoionization 
rate = 0.015 s~^, yielding (xhi) ~ 0.1. This value corre- 
sponds to the lower limit found by previous works (MH04; 
Wyithe & Loeb, 2004); in a forthcoming paper we will 
present results relative to other initial xhi values. In this 
work we do not account for spatial fiuctuations in the ul- 
traviolet background (UVB), which are believed to be in 

lower we can safely consider 30 per cent as an upper limit for the 
temperature increase due to helium photoheating. 
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place even in a highly ionized IGM at this epoch (Maselli 
& Ferrara 2005). The main implication of the assumption 
of a uniform background is a possible underestimate of the 
deviation from spherical symmetry of the physical Ifll re- 
gion. Accounting for a possible inhomogeneous background 
is expected to increase the dispersion of the physical radius 
of the HII region along different lines of sight. 

FiglU shows the xm distribution across the quasar loca- 
tion a.tt = tQ, the end of the RT simulation. The H ii region 
does not exhibit strong deviations from spherical symmetry. 
This result is not unexpected: the radiative energy density 
inside the H ii region during the early phases of the evo- 
lution is so large that clumps possibly responsible for flux 
anisotropics are completely ionized and made transparent. 
RT effects are instead apparent in the jagged ionization front 
(IF), causing the size of the H ii region to fluctuate along 
different LOS. We deflne the radius of the H ii region, Rd, 
along a given LOS as the distance from the quasar at which 
xhi > 10~^, marking the IF. The RT- induced scatter in the 
radius of the H ii region is seen in Fig. [2] via the probability 
distribution function (PDF) of Rd resulting from a sample of 
1000 LOS piercing the box through the quasar position. The 
mean value, (Rd) = 6.29±0.37 (1-cr), matches quite well the 
one derived from eq. [T] In addition, the uncertainty on Rd 
induced by RT effects is likely smaller than the experimental 
error on zq determinatior0. 

Next, we derived 1000 mock quasar absorption spec- 
tra along the same set of LOS used for Rd- The details of 
the adopted technique are given in Gallerani etal. (2006); 
in brief, each spectrum is characterized by a spectral reso- 
lution TZ — A/AA ~ 8000. To enable comparison with data 
each spectrum has been smoothed to 7?. = 4500 and Gaus- 
sian noise has been added, yielding a signal-to-noise ratio 
S/N = 50 (see e.g. Fan etal. 2006b). 

From these spectra we aim at deriving the observed 
H II region radius, Rf. In general, Rf Rd due to possi- 
ble effects of the Lya damping wing absorption arising from 
H I located outside the H ii region and to resonant absorp- 
tion from H I inside it. The deflnition oi Rf is somewhat 
arbitrary, as the transmissivity of the IGM at z ~ 6 is a mix- 
ture of dark gaps and transmission peaks (Fan etal. 2006a). 
As a consequence, the edge of the H ii region cannot be sim- 
ply identifled with the first point at which the transmitted 
flux drops to zero. Two different methods have been used so 
far in the literature: (i) Rf corresponds to the red side of the 
GP trougljfl, (ii) ^/ is identifled by the redshift at which the 
transmitted flux is > 0.1, when the spectrum is rebinned to 
AA=20A (Fan etal. 2006a). We have applied both meth- 
ods to derive Rf from our synthetic spectra and found only 
marginal discrepancies. The Rf PDF obtained from method 
(i) is shown in Fig. [2] (top panel) . From the Figure a large 



^ Accurate quasar redshift determinations are compromised by 
the systematic velocity offset between emission Unes of highly- 
ionized elements {i.e. C IV and Si IV ) and narrow lines probing 
directly the host galaxy (i.e. O III or CO molecular lines). Fan 
etal. (2006b) quantify an induced mean error in the measured 
redshift of Az ^ 0.02, corresponding to a proper distance of ~ 1.2 
Mpc at 2 ~ 6. 

^ As far as this work is concerned, we are taking into account the 
Lya GP trough. The analysis of the Ly/3 region could provide a 
larger Rf size. See Section 3 for further discussion. 
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Figure 2. Upper panel: Probability distribution function for Rd 
and Rf (physical units) using 1000 LOS through the simulation 
box. The offset between the two distributions quantifies the ap- 
parent shrinking (see text). Central panel: Illustrative template 
(dashed line) and absorbed (solid dark) spectra, along with the 
rijji density distribution (light gray) as a function of observed 
wavelength, for a representative LOS. Bottom panel: Contribu- 
tions to the total GP optical depth, t, from neutral hydrogen 
within (rr, dark) and outside [rd, light gray) the H II Region for 
the same LOS shown in the middle panel. 



offset between {Rd) = 6.29 Mpc and {Rf) = 4.25 Mpc is 
seen: i. e. the size of the H ii region extracted from the spec- 
tra is systematically underestimated. We refer to this effect 
as apparent shrinking. Also shown in Fig[2] (middle panel) 
are the template and absorbed spectra, along with the nni 
density distribution as a function of observed wavelength, 
for a representative LOS. 

The total GP optical depth r, responsible for the appar- 
ent shrinking, is the sum of two contributions: the damping 
wing absorption arising from H i outside the H ii regiory, 
and the resonant one, r^, from residual H i inside it. A de- 
tailed analysis of the mock spectra shows that, for xhi ~ 0.1, 
the Tr contribution to r is dominant. This can be appreci- 
ated from the lower panel of Fig[5] where Tr and Td are 



* If the quasar H II region is embedded in a partially neutral 
IGM, the H I outside the H II region, can produce significant ab- 
sorption at wavelengths that in physical space correspond to the 
ionized region (see e.g. Madau & Rees, 2000). 
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Figure 3. Likelihood contours (68 and 95 per eent confident 
limit) for the Rf distribution in the (xHijig) plane. The cross 
and triangle indicate the most likely values {xiii,tQ) for Rf and 
Rd distributions, respectively. The square is the same quantity 
derived from the observed quasar sample. 

plotted separately along the same representative LOSQ The 
substantial contribution of resonant absorption results from 
the increase of the average xm with physical distance from 
the quasar due to flux geometrical dilution and attenuation. 
Close to the edge we flnd Tr ~ 400, on average. 

The apparent shrinking introduces a mean systematic 
underestimate of the physical H ii region size, Rd, by AR — 
{{Rd - Rf)/Rd) = 0.32. Note that the amplitude of AJ? is 
well above errors induced by RT effects and uncertainties 
in the quasar parameters. Ai? has a considerable dispersion 
around the mean value above, mostly due to the large fluc- 
tuation of Rf along different LOS (see Fig. [S]). In addition, 
there is no speciflc correlation between Rd and Rf along 
different LOS. Both these effects can be understood from 
the fact that Rf depends on the Lya optical depth, which 
in turn is much more sensitive than the ionizing continuum 
opacity to tiny fluctuations of xm inside the bubble. In or- 
der to properly predict such fluctuations, is very important 
to perform accurate RT calculations. 

In practice, though, the lack of correlation between Rd 
and Rf makes it very difBcult to precisely derive Rd from 
the observed spectrum, making it crucial to quantify the re- 
liability of Xm constraints obtained from Rf measurements 
in high-2 quasar spectra. To address this issue, we have per- 
formed a maximum likelihood analysis of our results. Given 
the Rf distribution derived from our mock spectra by ap- 
plying method (ii), we have calculated its likelihood func- 
tion (LF) to match a Gaussian distribution, Q{R), of mean 
value given by eq.[T] and r.m.s. equal to that of Rf. We im- 

® As Td is oc xhi outside the HII region, it only becomes com- 
parable to Tr when xhi 0.5, and overcomes in the case of a 
complete neutral IGM. 



pose the following priors: iV^- — 5.2 x 10^*^ s~^, xni £ [0, 1], 
tq G [10^, 10*] yr. For the quasar lifetime tq we assume an 
a lognormal prior with mean equal to tq = lO'^ yr; uncer- 
tainties in N-^ can be absorbed into tq variations. We have 
rejected (xHijig) pairs yielding Rd values outside the extent 
of the Rf distribution. The LF maximum, {xiii,tq), identi- 
fies the most likely values inferred from a given sample of 
observed quasar spectra. We find {xm^tq) = (0.34, lO'^yr), 
point marked by a cross in FigO however, these values 
lay 2-(j away from the actual values used in the simulation 
{x-m,tq) ~ (0.1, lO^yr), i.e. those that we would like to re- 
cover. As a sanity check we have repeated our maximum 
likelihood analysis on the Rd distribution, and indeed we 
find that the maximum coincides with the simulation values 
(triangle in Fig. |3J . We conclude that the apparent shrink- 
ing effect induces an overestimate of the xbj by a factoi[f| 
~ 3. It is worth noting that, although our Rf distribution 
has been drawn from a simulation of a single quasar with 
fixed intrinsic properties {N^, tq, zq), the range of accept- 
able xm values is quite large. For example, fixing tq to the 
simulation value, 10^ yr, still allows xbj > 0.1,0.07 to a 1 - 
a,2-a confidence level. Possible dispersion in the intrinsic 
properties of the quasar sample are likely to make the xm 
determination even more difficult. 

We have also applied the maximum likelihood analy- 
sis to a sample of 6 observed QSOs, having zq G [6.0,6.2], 
whose spectra have been studied by Fan etal. (2006a). Such 
authors give the measured radii scaled to a reference com- 
mon quasar absolute magnitude M1450 = —27. We have then 
calculated the LF of the sample of observed spectra to match 
a Gaussian with mean value given by eq[T]where A^.^ is scaled 
to[3| to M1450 = —27, retaining the observed sample lumi- 
nosity dispersion. 

In this case the LF maximum is {xm,tq) — (0.06, 10^), 
shown as a square in Fig[3] By taking into account the over- 
estimate of this result due to the apparent shrinking effect, 
and the uncertainties on the location of the maximum in- 
duced by the wide range of acceptable xm values, we find 
that our study slightly favors a mostly ionized universe at 
z ~ 6.1. By varying Nj in the range [2.7, 7.7] x 10^® suggested 
by MH04, and scaling the theoretical radii accordingly, our 
results are always consistent with xm < 0.2. 



3 SUMMARY AND DISCUSSION 

In this Letter we have discussed the robustness of con- 
straints on the IGM H i fraction inferred from the extent of 
H II regions around high redshift quasars by means of their 
absorption spectra. 

We have performed a combination of state-of-art mul- 
tiphase SPH and 3D RT simulations to reliably predict the 
properties of a typical high-z quasars H 11 region {e.g. ex- 
tent, geometrical shape, inner opacity). We have found that 
RT effects do not induce strong deviations from spherical 
symmetry. The RT-induced dispersion in the H 11 region size 

^ This factor could depend slightly on xjji and, to a lesser extent 
on Nj. 

We use the scaling: N.y = N^fi X io(-27-Mi45o,o)/2.5^ -yVg 
A.y,0 = (5.2 ±2.5) X 105*^8-1 and Mi45o,o = -27.2 which are the 
estimated values for the z = 6.28 QSO SDSS J1030+0524. 
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along different LOS is in fact of the order of roughly 6 per 
cent of the mean radius which is likely smaller than the typ- 
ical error induced on Rd estimates by uncertainties in the 
quasar redshift determinations. 

By deriving and analyzing mock spectra through the 
simulated quEisar environment we have found that the 
H II region size deduced from quasar spectra, Rf, typically 
underestimates the physical one by 30 per cent. 

The fact that the observed H ii region sizes can substan- 
tially underestimate the size of the region impacted by the 
ionizing radiation of the quasar was already noted by Bolton 
& Haehnelt (2006), but our results gives the first quantita- 
tive estimate of the above underestimation, to which we refer 
to as apparent shrinking. It is worthwhile to notice here that 
the resolution of our simulations is not sufficient to properly 
resolve the IGM clumping which could result in an under- 
estimate of RT effects. To this respect our results can be 
considered a lower limit for the amplitude of the apparent 
sfirinking effect: resolving the high density clumps would re- 
sult in fact in a higher inner opacity which is s expected to 
be the origin of a stronger apparent shrinking. 

The amplitude of the apparent shrinking effect could be 
decreased by defining Rf as the red side of the GP trough 
in the Ly/3 region (ii^). A clear offset between the red side 
of the GP trough in Lya and Ly/3 has been observed in the 
spectrum of the QSO SDSS J1030-I-0524 (White etal. 2003, 
MH04. Recently, Bolton & Haehnelt (2006) have argued that 
the ratio R^/Rf has a well-defined trend with xhi, which 
could be exploited to constrain such quantity. However, the 
robustness of this result is uncertain given the small number 
of LOS used in the study and further investigation is needed. 

In addition our analysis shows that, in the case of 
xm = 0.1, the apparent shrinking is almost completely due 
to resonant absorption of residual H i inside the ionized bub- 
ble. This contribution is highly fluctuating along different 
LOS, resulting in a large dispersion of the observed radii 
distribution, which is mainly due to the fact that Rf de- 
pends on the inner Lya opacity which is highly sensitive to 
riHi fluctuations. This result implies that in order to prop- 
erly account for these fluctuations, detailed RT calculations 
arc required; the present work is unique in this respect, be- 
ing the only study which accounts for the detailed effects of 
3D RT. 

The maximum likelihood analysis we have performed 
on a sample of 1000 mock spectra shows that the apparent 
shrinking effect induces an overestimate of the xhi by a fac- 
tor w 3, if the IGM is only partially ionized {xm = 0.1). 
Moreover, by applying the above analysis to a sample of ob- 
served QSOs, we conclude that our study favors a mostly 
ionized universe at 2: ~ 6.1 (xhi ^ 0.06), a somewhat differ- 
ent conclusion with respect to the previous results given in 
MH04 and Wyithe & Loeb (2004) and Wyithe etal. (2005), 
which found a lower limit for xm of 0.2 and 0.1 respec- 
tively. Nevertheless, Bolton & Haenhelt (2006) posted as an 
outcome of their analysis that the sizes of the H 11 regions 
in the higher redshift quasar spectra are consistent with 
a significantly neutral surrounding IGM, as well as with a 
highly ionized IGM. Furthermore, as mentioned in the In- 
troduction, independent studies by Yu & Lu (2005) and Fan 
etal. (2006a) give evidence for much lower values of xm- 
It is worth noting that our measurement agrees with the 
independent determination by Fan etal. (2006) who found 



xui — 1.3 X 10^'^, based on an analysis of the HIT region size 
evolution. 

Uncertainties remain due to the fact that the range of 
acceptable xm values associated to a given sample of mea- 
sured radii is quite large, e.g. our Rf sample, having as most 
probable value xhi = 0.34, still allows xm > 0.1, 0.07 to a 
l-a,2-a confidence level. This suggests that measurements of 
the H II size in quasar spectra can only provide rough con- 
straints on Xm, as long as the knowledge of intrinsic prop- 
erties of observed QSOs remains incomplete. 

It is worth noting that quasar spectra could contain ad- 
ditional useful information on xm- MH04 constrained xm 
from an analysis of the SDSS ,11030+0524 absorption spec- 
trum, which is completely independent from constraing xhi 
using the radius extent. In this case the xhi determination 
is based on the difference between the Lya and Ly/3 optical 
depths and it is free from the uncertainties related to the 
apparent shrinking effect. It is interesting to notice that, dif- 
ferently from our results, MH04 find that the slightly smaller 
size of the H 11 region extent inferred from the Lya and Ly/3 
opacity is due to the presence of a damping wing absorption 
produced by the H i in the IGM unperturbed by the quasar. 
This difference can bo duo to the fact that MH04 compute 
the neutral hydrogen fraction at the equilibrium with the 
quasar ionizing flux which is diluted geometrically, without 
taking into account the further attenuation due to the inner 
opacity. As a consequence the inner opacity itself is under- 
estimated. We found instead that the inner opacity, whose 
correct determination requires detailed 3D RT calculations, 
must be taken into account for a reliable determination of 
the apparent shrinking effect. A significant underestimate of 
the inner (resonant) opacity would require a higher xm out- 
side the H II region in order to reproduce the observed dif- 
ference between the apparent shrinking in the Lya and Ly/3 
flux. In any case, though, it is necessary to await for a larger 
sample of observed quasar and simulate a larger parameters 
space (xm, tq, N.-,) before being able to directly compare 
results coming from different approaches in the modelling as 
well as in the statistical analysis technique adopted. 
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